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Landslide inventories show that the statistical distribu- 
tion of the area of recorded events is well described by a 
power law over a range of decades. To understand these 
distributions, we consider a cellular automaton to model a 
time and position dependent factor of safety. The model is 
able to reproduce the complex structure of landslide distri- 
bution, as experimentally reported. In particular, we inves- 
tigate the role of the rate of change of the system dynamical 
variables, induced by an external drive, on landslide mod- 
eling and its implications on hazard assessment. As the 
rate is increased, the model has a crossover from a critical 
regime with power-laws to non power-law behaviors. We 
suggest that the detection of patterns of correlated domains 
in monitored regions can be crucial to identify the response 
of the system to perturbations, i.e., for hazard assessment. 

Introduction As for earthquakes and forest fires, there 
is a compelling evidence that the landslide frequency-size 
distributions are power- law functions of the area [Turcotte 
et al., 2002]. The presence of these broad distributions has 
crucial consequences on both the basic understanding of 
these phenomena and the practical and relevant purposes, 
such as the evaluation of natural hazards. Here, we intro- 
duce a cellular automaton that is aimed at modeling the 
general features of landslides, and is focused on the dynam- 
ical evolution of a space and time dependent factor of safety 
field. This model is very simple, but it is able to give a 
comprehensive picture of the avalanching phenomena and 
to reproduce some well-known properties of landslide distri- 
butions. 

Several authors invoked the paradigm of self-organized 
criticality (SOC) [Bak et ai, 1987; Jensen, 1998; Turcotte, 
1999] to explain landslide distributions [Turcotte et al, 2002; 
Pelletier et al., 1997; Hertgarten and Neugebauer, 2000]. Al- 
though the "critical" nature of the present phenomenon is 
not yet assessed and many authors believe that deviations 
from power-law appear to be systematic for small land- 
slides data [Stark and Hovius, 2001; Pelletier et al., 1997; 
Brardinoni and Church, 2004; Malamud et al., 2004], several 
regional landslide inventories records show robust power- 
law distributions of large events with an exponent around 
a ~ 2.5 [Turcotte et al., 2002], ranging approximately from 
a ~ 1.75 for rockfalls to a ~ 2.8 for mixed landslides (see 
[Dussauge et al., 2003; Faillettaz et al., 2004] and refer- 
ences therein). The universality of such an exponent is still 
debated (see [Turcotte et al., 2002; Dussauge et al., 2003; 
Faillettaz et al., 2004; Malamud et al., 2004] and references 
above), and its reported values are far from the one in the 
original "sandpile model" [Bak et al., 1987], where a ~ 1.0. 
Recently, the reported values of a have been obtained by 
introducing two-thresholds mechanisms in models that re- 
late landslide dynamics both to SOC [Pelletier et al., 1997; 
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Hertgarten and Neugebauer, 2000] and to non SOC cellu- 
lar automata [Faillettaz et al., 2004]. Furthermore, models 
based on F [Malamud et ai, 2004] or Pareto [Stark and Hov- 
ius, 2001] distributions have been proposed. 

Here, we consider a model inspired to an anisotropic ver- 
sion of the Olami-Feder-Christensen (OFC) [Olami et al, 
1992] cellular automaton and subject to a finite driving rate 
[Hamon et al, 2002]. The model describes the evolution of 
a space and time dependent factor of safety, which is in- 
vestigated for the first time in the present framework. In 
particular, we outline the essential role played by the rate 
of change of the system dynamical variables (variation of 
pore water pressure, lithostatic stress, cohesion coefficients, 
etc., [Helley et al., 2004; Iverson et al, 2000]) induced by 
external triggers. We find the model to be at the edge of 
the SOC limit. Actually, such a limit, which is achieved 
only in the asymptotic condition of vanishing driving rate, 
is hardly attainable in a real landslide process. The model 
is able to reproduce power-law distributions with exponents 
very close to the observed values. Power-laws are robust 
even though their exponent smoothly depends on system 
parameters (e.g., time derivative of the factor of safety and 
its dissipation level, see below). In this sense, although the 
SOC paradigm to some extent may be applied to landslides 
[Turcotte et ai, 2002], the idea of universality, within this 
model, must be restricted to the shape of the frequency-size 
distribution rather than to its exponent, as can be deduced 
from some catalogues [Dussauge et al, 2003; Faillettaz et al, 
2004]. Finally, in presence of strong driving rates we find 
that the model has Gaussian behaviors. We examine below 
the implications of our results on hazard assessment. 

The Model The empirical Mohr-Coulomb failure crite- 
rion establishes that landslides occur when the shear stress 
exceeds a maximum value, which is given by T max — c + 
(a — u)tan<j!>, with a the total normal stress, u the pore- 
fluid pressure, <j> the angle of internal friction of the soil 
and c the cohesional (non-frictional) component of the soil 
strength [Terzaghi, 1962]. In literature, the factor of safety, 
FS, against slip is defined by the ratio of the maximum 
shear strength r max to the disturbing shear stress r 

T 

If FS > 1 resisting forces exceed driving forces and the slope 
remains stable. Slope failure starts when FS = 1. Although 
the practical determination of FS is difficult, simple one- 
dimensional infinite-slope models can quantify how c, u and 
<f> influence the Coulomb failure and show that the ground- 
water term has the most widely ranging influence [Iverson 
et ai, 1997]. Traditional models generally treat soils and 
rocks as continuous porous media that obey to the Darcy's 
law. Actually, field evidence indicates that the hydrology 
of some natural slopes is strongly influenced by discontinu- 
ities such as fractures and macropores. In practice, obser- 
vations of large spatial and temporal fluctuations of water 
flow, within slopes at different sites, support the assertion 
that water-flow paths and permeability continually change 
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within the slopes and also during the failure, providing dif- 
ferent local values of the pore pressures and of the cohesion 
[Iverson et al, 1997; Helley et al, 2004; Iverson et al, 2000]. 

In order to take into account the complex non- 
homogeneous structure of a slope in the above failure condi- 
tion, we consider a site and time dependent factor of safety, 
FS. In particular, we approximate a natural slope by a two- 
dimensional (square) grid and define on each cell, i, of the 
lattice a local variable d — 1/FSi. Such a local inverse fac- 
tor of safety is the fundamental dynamical variable of our 
model. The presence of diffusion, dissipative and driving 
mechanisms acting in the soil, such as those on the water 
content, inspires the dynamics of our model, which is de- 
fined by the following operations. Starting from a random 
and "stable" initial configuration (e^ < 1 Vi), the system is 
subject to changes caused by some external trigger, as for 
instance a uniform rainfall, and the values of e; on each cell 
of our grid change at a given rate v, — > e, + v. For the 
sake of simplicity, we consider here only a uniform driving 
rate, v, but different choices can be made to simulate the 
effect of different hydrologic and external triggering mecha- 
nisms. The model is driven as long as e, < 1 on all sites i. 
Then, when the generic site i becomes unstable (i.e., over- 
passes the threshold, a > 1), it relaxes with its neighbors 
according to the rule: 

6z > 0; 6nra > e nn ~\~ fnn&i> 

where the index nn denotes the nearest neighbors of site i 
and /„„ is the fraction of e; toppling on nn (after failure 
we set 1/FS = for simplicity, as any other finite level 
would work [Jensen, 1998]). This kind of chain relaxations 
("avalanches") is considered to be instantaneous compared 
to the time scale of the overall drive and it lasts until all sites 
are below threshold. The model is said to be conservative 
if C — X/nn f nn ~ Since many complex dissipative phe- 
nomena (such as evaporation mechanism, volume contrac- 
tions, etc. [Fredlund and Rahardjo, 1993]) contribute to a 
dissipative stress transfer, we consider the non-conservative 
case C < 1, which is different from previous landslide mod- 
els [Pelletier et al, 1997; Hertgarten and Neugebauer, 2000]. 
Since gravity individuates a privileged direction, we consider 
an anisotropic model where the fraction of e» moving from 
the site i to its "downward" (resp. "upward") neighbor on 
the square lattice is fd (resp. /„), as fi = f r is the fraction 
to each of its "left" and "right" neighbors. In particular, we 
assume f u < fd and fi = f r < fd- This choice of parameters 
is made in the attempt to sketch the complex relaxation pro- 
cesses occurring in a slope failure. The conservation level, 
C, and the anisotropy factors, f's, which we assume to be 
uniform, are actually related to the local soil properties (e.g., 
lithostatic, frictional and cohesional properties), as well as 
to the local geometry of the slope (e.g., its morphology). 
The rate of change of the inverse factor of safety, v, which is 
induced by the external drive (e.g., rainfall) and is related 
to soil and slope properties, quantifies how the triggering 
mechanisms affect the time derivative of the FS field. 

Numerical Results We consider a 64 x 64 square lat- 
tice, implementing both cylindrical (open along the verti- 
cal axis and periodical along the horizontal axis) and open 
boundary conditions, which do not give appreciable differ- 
ences. Once the system has attained a stationary state in 
its dynamics, we study the probability distribution, P(s), of 
avalanches of size s. During a run (we treat statistics of over 
10 9 events per run) the conservation level C and the rate, v, 
are kept fixed. Examples of the frequency-size distribution 
of avalanches, P(s), are shown in figures 1 and 2. In figure 
1, the different curves correspond to different values of the 
rate v, for C = 0.4, and in figure 2 to different values of C, 
for v = 5 ■ 10 -3 . 

In the limit of very small driving rate, i.e., v — > 0, the 
distribution of events, P(s), exhibits the typical SOC struc- 
ture (see figure 1): a power law characterized by a critical 



exponent a, P(s) ~ s~ a , in agreement with the experi- 
mental evidence on medium and large landslides [Turcotte 
et al, 2002; Pelletier et al., 1997; Hertgarten and Neuge- 
bauer, 2000; Faillettaz et al, 2004; Brardinoni and Church, 
2004; Malamud et al, 2004; Dussauge et al, 2003], followed 
by a size dependent exponential cutoff [Jensen, 1998]. By 
increasing the rate v, the power-law regime shifts towards 
larger sizes and at some point the probability distribution 
apparently shows a maximum for a value s* . There are 
two regimes: for large landslides (s > s*) the above struc- 
ture P(s) ~ s~ a is found, while for small events (s < s*) 
an increasing function of s is observed. Such a complex 
structure is absent in SOC models, instead a maximum is 
found in landslide inventory maps for small landslide data, 
although there is no consensus about the nature of such a 
feature [Stark and Hovius, 2001; Brardinoni and Church, 
2004; Malamud et al, 2004]. The values of the power-law 
exponent, a, by varying the rate, v, and C are very close 
to those experimentally found [Turcotte et al, 2002; Dus- 
sauge et al, 2003]. As in the original isotropic OFC model 
[Olami et al, 1992], the critical exponent decreases with the 
level of conservation, C (see inset in figure 2). The value 
of a slightly changes with the anisotropic ratios fd/ fd and 
fu/ fd, except when they get too small [Piegari et al, 2005] 
where, as found also in other models [Amitrano et al, 1999; 
Amitrano, 2003; Faillettaz et al, 2004], the event size dis- 
tribution is considerably modified. The power-law regime is 
crucially robust to changes in system parameters. For in- 
stance, in the case of figure 1 it can be found for v up to 
approximately 10 -2 , all over the range C G [0.4,0.8]. It is 
worth noting that the a values here obtained are comparable 
to those found in models of failure in fiber bundle [Hemrner 
et al, 1992; Hansen et al, 1994; Hidalgo et al, 2002]. 

As it can be seen in figure 1, a further increase of the driv- 
ing rate (above v ~ 10 -2 ) causes a crossover to a markedly 
different regime where power-laws are no longer apparent 
and a bell shaped (Gaussian) distribution emerges, whose 
peak shifts towards larger sizes and shrinks up. Such a be- 
havior is to be expected since for strong driving rates all 
internal correlations are washed out. 

Summarizing, the conservation level, C, and the time 
derivative of 1/ FS, v, turn out to be important to determine 
landslide probability distributions: in the limit v — ► 0, the 
model is indeed in the SOC class; for small but finite v the 
system is at the edge of SOC and the critical behaviors are 
still largely observed; finally, as v gets large enough, Gaus- 
sian properties are found. Thus, in the small v regime, our 
model reproduces the general properties of existing catalogs 
and can help interpreting them, while in the large v regime 
it foresees a different class of behavior. Nevertheless, for the 
sake of clarity, we have considered the simple case where the 
rates v and C are uniform. Thus, the distributions of figure 
1 may be not directly comparable to landslide inventories, 
which gather events with non-uniform driving rates. 

Pictures of a typical "avalanche" in the different regimes 
discussed above are plotted in figure 3 (upper panels) with 
the corresponding values of the factor of safety, FSi, on the 
model grid after the avalanche (lower panels). The snap- 
shots in the upper panels, taken in two systems driven at 
different rates (left v = 5 ■ 10~ 3 , right v = 5 ■ 10~ 2 ), show 
a typical event with size s — 230 (such a value is chosen 
because it has approximately the same probability in the 
two cases, see figure 1): the system on the left is in the 
power-law regime; the one on the right is in the non power- 
law regime. The difference of the avalanche geometry in the 
two cases is impressive. Domino effects are crucial to de- 
termine a "catastrophic" event when the system is governed 
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by a power-law statistics, where a huge compact landslide 
is present [Pietronero and Schneider, 1991] with a typical 
size of the order of the system size (left-upper panel). Con- 
versely, large events are expected at higher v (where indeed 
the average size (s) is much larger than 230), but in such a 
regime a typical event with s = 230 is made of many tiny 
unconnected avalanches (summing up to s — 230). 

Interestingly, even though the P(s) is very different in the 
two cases, the probability distribution, P(FS), of the spa- 
tial values of FS on the grid has a similar Gaussian shape, 
with comparable averages (FS) and fluctuations (AFS 2 ) 
({FS) = 2.20 and (AFS ) = 0.16 for v = 5 • 10~ 3 ; 
(FS) = 2.53 and (AFS 2 ) = 0.06 for v = 5 • 10~ 2 ) lay- 
ing far above the instability threshold FS = 1. Thus, a 
measure of just an average safety factor on the investigated 
area could provide only very partial information about the 
statistics governing landslide events. 

The origin of the striking difference of the P(s) in the 
two considered cases traces back to the relative extension 
of spatial correlations of the factor of safety, FS, which is 
derived from the correlation function: 

r( ~ = (FS(r)FS(r + z)) - (FS(r)) 2 
( ' (FS(r) 2 ) - (FS(f)> 2 V ' 

where FS(x) is FS at position x (here we take z along 
the direction of the slope) and the average is over the sys- 
tem sites. As it is well-known [Jensen, 1998], we find that 
C(z) oc cxp(— «/£), where £ is the spatial correlation length 
of FS. The value of FSi on site i is shown in the lower pan- 
els of figure 3, in gray scale, for the same cases pictured in 
the upper panels: patterns of large correlated areas (regions 
with similar values of FSi, i.e., the same color) are appar- 
ent in the left bottom panel and, in practice, absent in the 
right one. In the power-law regime (e.g., v = 5 • 10 -3 ), the 
correlation length, £, is of the order of the system size; thus, 
even a very small perturbation (say, a drop of water) at one 
single point can trigger huge system responses. Instead, in 
the non power law regime (e.g., v = 5- 10~ 2 ) large-scale cor- 
relations are absent; here, large events trivially occur just 
because the strong external driving rate makes likely that 
many cells simultaneously approach the instability thresh- 
old. The detection of patterns of correlated domains (i.e., 
the size of £) in investigated areas results, thus, to be a 
crucial tool to identify the response of the system to pertur- 
bations, i.e., for hazard assessment. 

Conclusions To summarize, we have investigated a 
continuously driven anisotropic cellular automaton model 
for the characterization of landslide size distributions. The 
model may help in interpreting the general behaviors ob- 
served in real systems. In particular, we have found that 
different values of the driving rate give rise to different sta- 
tistical distributions of events. The determination of cor- 
related domains in the factor of safety becomes crucial for 
landslide classification and, consequently, for hazard assess- 
ment. 
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Figure 1. The probability density distribution, P(s), of 
avalanches of size s is plotted for the shown values of the 
time derivative of the inverse factor of safety, v (model 
size L 2 — 64 x 64, conservation level C = 0.4, anisotropy 
coefficients f u /fd = 2/3 and fi/fd = 5/6). The power 
law P(s) ~ s~ a found in the limit v — ► is partially pre- 
served by increasing v up to a point where a bell shaped 
behavior is clearly observed. 
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Figure 2. The probability distribution, P(s), of 
avalanches of size s is plotted for the shown values of 
level of conservation, C (y = 5 • 1CU 3 , other param.s as in 
figure 1). The inset shows the exponent a (confidence 
interval 95%) of the power-law fit as a function of C. 



Figure 3. Top panels: Two snapshots of a typical 
landslide event of size s = 230 on our 64 x 64 grid, in 
two cases with a driving rate v = 5 • 1CU 3 (left figure) 
and v = 5 • 1CU 2 (right figure). The cells marked in black 
are those which reached the instability threshold. Bot- 
tom panels: The pictures plot the local value of the 
factor of safety, FS, corresponding to the stable config- 
urations reached after the avalanches shown in the up- 
per panels. The FS values have been associated to ten 
levels of color from white to black, in order to measure 
the distance of a cell from its instability condition: the 
darker the color, the farther is the cell from the instability 
threshold. In the panels it is possible to recognize as dark 
areas the avalanches shown in the corresponding upper 
grids. In particular, dark areas are related to previous 
landslide events, as the lighter areas indicate regions of 
future events. In the left figure large correlated regions 
(compact areas with same color) are observed, whereas 
their size is small in the right figure. 



